Microscopic Model and Phase Diagrams of the Multiferroic Perovskite Manganites 
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Orthorhombically distorted perovskite manganites, RMxiO^ with R being a trivalent rare-earth 
ion, exhibit a variety of magnetic and electric phases including multiferroic (i.e. concurrently mag- 
netic and ferroelectric) phases and fascinating magnetoelectric phenomena. We theoretically study 
the phase diagram of .RMnOs by constructing a microscopic spin model, which includes not only the 
superexchange interaction but also the single-ion anisotropy (SIA) and the Dzyaloshinsky-Moriya 
interaction (DMI). Analysis of this model using the Monte-Carlo method reproduces the experimen- 
tal phase diagrams as functions of the 7?-ion radius, which contain two different multiferroic states, 
i.e. the afc-plane spin cycloid with ferroelectric polarization P\\a and the fee-plane spin cycloid with 
P\\c. The orthorhombic lattice distortion or the second-neighbor spin exchanges enhanced by this 
distortion exquisitely controls the keen competition between these two phases through tuning the 
SIA and DMI energies. This leads to a lattice-distortion-induced reorientation of P from a to c in 
agreement with the experiments. We also discuss spin structures in the A-type antiferromagnetic 
state, those in the cycloidal spin states, origin and nature of the sinusoidal coUinear spin state, and 
many other issues. 

PACS numbers: 75.80.-|-q, 77.80.-e, 75.30.Gw, 75.47.Lx 



I. INTRODUCTION 



Rare-earth manganites with orthorhombically dis- 
torted perovskite structure, i?Mn03 with R being a triva- 
lent rare-earth ion, have been subject to intensive stud- 
ies since the multiferroic phases, in which magnetism 
and ferroelectricity simultaneously emerge, were found 
in some of these materials Jj, i2(] . This class of materials 
exhibits a variety of magnetic and electric phases as a 
function of ionic radius of the R ion. The i?-site varia- 
tion controls magnitude of the orthorhombic lattice dis- 
tortion (the GdFeOa-type distortion) see Fig. [T]^a) 

, i.e. tilting angles of the MnOg octahedra become 

larger and the perovskite lattice is more significantly dis- 
torted with a smaller-sized R ion. An experimentally 
obtained phase diagram in Refs. 0, 13 exhibits A-type 
antiferromagnetic [AFM(A)] ground states in the weakly 
distorted materials with _R=La, Pr,..., Eu, Gd, while E- 
type antiferromagnetic [AFM(E)] ground states in the 

strongly distorted materials with i?=Ho,..., Yb, Lu 

see Fig. [Hb). Sandwiched by these two regions, a spiral 
spin order with concomitant ferroelectricity is observed in 
the moderately distorted materials with i?=Tb and Dy. 
Anticipated coupling between spins and electric dipoles 
(magnetoelectric coupling) makes these materials inter- 
esting because of possible technical applications. 

In TbMnOs, for example, a sinusoidal coUinear order 
of Mn spins occurs at the Neel temperature Tj^^" ~ 41 
K d, @| • In this phase, the Mn spins are aligned along the 
h axis with an incommensurate propagation wave vector 
9m" =(0, 0.28, 1) in the Pbnm orthorhombic unit cell. Be- 
low Tc ^ 28 K, ferroelectricity shows up along the c axis 
concomitantly with a magnetic transition into a trans- 
verse spiral (cycloidal) spin order with Mn spins rotating 



within the be plane 7] • Upon further decreasing temper- 
ature, ordering of the /-electron moments on the rare- 
earth (Tb) ions takes place at '--^ 7 K with a differ- 
ent propagation wave vector ~(0, 0.42, 1). DyMnOa 
shows similar orderings and transitions with Tj^" ~39 K, 
qT ^(0, 0.36, 1), Tc -19 K, rj« ~ 5 K, and q« ~(0, 
0.5, 1) 8]. 

The ferroelectricity in these materials can be micro- 
scopically explained by the spin-current model [in other 
words, inverse Dzyaloshinsky-Moriya (DM) mechanism] 
proposed by Katsura, Nagaosa and Balatsky (KNB) Q. 
Overlaps of electronic wave functions between adjacent 
two atomic sites (i and i -I- 1) with mutually canted spins 
{Si and Si+i) can generate a local electric polarization. 



p,; = Ae,^,+i X (5j X St+i), 



(1) 



where Bij^i denotes the unit vector connecting these 
two sites, and A is a constant determined by the spin- 
exchange and spin-orbit interactions. The cycloidal spin 
order in i?Mn03 is expected to induce a uniform spon- 
taneous polarization P as a sum of the local polar- 
izations Pi in the direction perpendicular to the spi- 
ral propagation vector (||6) and the spin-helicity vector, 
hs = J2i'^i ^ '^i+i- Consequently, the ferroelectric po- 
larization parallel to the a axis (Pa) is expected in the 
a6-plane cycloidal spin structure, while in the 6c-plane 
cycloidal spin structure, that parallel to the c axis (Pc) is 

expected see Fig. [2j Direct evidence for the inverse 

DM mechanism was recently demonstrated by controlling 
the spin-helicity vector in TbMnOa with an external elec- 
tric field [10]. Similar theoretical results have also been 
obtained independently in Refs. [U[l3l- In particular, 
Sergienko and Dagotto discussed that a spiral spin or- 
der can induce uniform shifts of the oxygen ions through 
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AFM(A)+WFM AFM(E) 

FIG. 1: (Color online) (a) Perovskite structure with GdFeOs-type distortion. The unit cell contains four Mn ions, which 
are referred to as Mn A, Mn B, Mn C, and Mn D, respectively, (b) Experimentally obtained magnetoelectric phase diagram 
of -RMnOs in plane of temperature and ionic ii-site radius reproduced from Refs. 0, [31 . The region between Gd and Tb 
was studied using a solid solution system Gdi-2,Tba;Mn03 [1^. Averaged 7i-site radii for the solid solutions are deduced 
by interpolation. Insets show spin configurations of the A-type and E-type antiferromagnetic [AFM(A) and AFM(E)] states. 
Along the c axis, spins stack antiferromagnetically. PE and FE denote paraelectric and ferroelectric phases, respectively. In the 
AFM(E) phase, possible ferroelectricity due to an exchange-striction mechanism was theoretically proposed [2ll. [2^ [23l|. and it 
was confirmed experimentally 24, 25, 26]. (c) Experimentally obtained magnetoelectric phase diagram of a solid solution system 
Eui-2,Ya;Mn03 in plane of temperature and Y concentration x reproduced from Ref. Insets show spin configurations of the 
a6-cycloidal, 6c-cycloidal, and sinusoidal collinear states. Spin configuration of the AFM(A) state with weak ferromagnetism 
(WFM) seen along the a axis is also shown. 



the spin-lattice coupling via the DM interactions, lead- 
ing to the electric polarizations given by Eq. ([T]). Recent 
first-principles calculations revealed that the Sergienko- 
Dagotto type lattice-mediated contribution is dominant 
for the ferroelectric polarization relative to the KNB type 
pure-electronic contibution [3, US • 

While only the Pc phase with 6c-cycloidal spin struc- 
ture is observed in the single-rare-earth compounds like 
TbMnOa and DyMnOa, the Pa phase with a6-cycloidal 
spin structure can be observed on the verge of the 
AFM(A) phase when we continuously control magnitude 
of the GdFeOa-type distortion using a solid solution sys- 



tem Gdi_a;Tba;Mn03 [T^. Moreover, reorientation of the 
polarization occurs from Pa to Pc with increasing Tb con- 
centration X see Fig. (Ub). This phenomenon is as- 
cribed to a spin-cycloidal-plane flop from the ab plane to 
the be plane [l6j . 

The reorientation of the electric polarization is 
also observed in another solid solution system 
Eui-ajYjjMnOs [13, [Hi. Phase diagrams of this 
system were experimentally studied for compositions 
^ 2; < 0.5 as functions of temperature and Y concen- 
tration x [13, [H, [H, [131 see Fig. ^c). According 

to these experiments, the ground state changes ap- 
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FIG. 2: (Color online) (a) Relationship between the spin- 
helicity vector hs = X^i ^ •S'l+i and the spontaneous elec- 
tric polarization P in the ab-cycloidal spin structure expected 
horn the spin-current model [Til . and (b) that in the 
6c-cycloidal spin structure. 



proximately at x ~ 0.2-0.3 from the canted AFM(A) 
state without long-range ferroelectric order towards the 
presumably a6-cycloidal spin state with Pa ■ For higher Y 
concentrations of a; > 0.4, the orientation of polarization 
spontaneously changes from Pc at higher temperatures 
towards Pa at lower temperatures. In addition, the 
regime of Pc phase increases with increasing x, which 
resembles the polarization flop with increasing x in the 
Gdi_a;Tba;Mn03 system. These results indicate that 
the polarization flop occurs both thermally and by the 
increase of GdPeOa-type distortion. These phenomena 
are caused purely by interactions among the Mn Sd 
spins since the Eui_a;Yj;Mn03 system is free from the 
influence of /-electron moments on the i? ions. 

The above experimental results tell us that the 
manganite system exhibits the AFM(A), a6-cycloidal, 
6c-cycloidal, and AFM(E) phases successively as the 
GdFeOs-type distortion increases. In particular, the 
cycloidal-plane flop between ah and he is quite important 
since the plane is related to the orientation of electric 
polarization. According to Ref. 3, enhanced tilting of 
the MnOg octahedra leads to an increase of the second- 
neighbor antiferromagnetic exchanges (J2), which com- 
pete with the nearest-neighbor ferromagnetic exchanges 
(Ji) in the ah plane. Resulting weakening of the effec- 
tive in-plane ferromagnetic interactions causes reduction 
of the Neel temperature in the AFM(A) phase as the 
(averaged) i?-site radius decreases. We can attribute the 
emergence of the long-wave-length antiferromagnetic or- 
ders or that of the spiral spin orders to magnetic frustra- 
tion caused by the above competition. The spiral spin 
order can be reproduced by a simple J1-J2 Heisenberg 
model [27I. [28I, [29} . Within this model, however, there 
exists degeneracy and the cycloidal plane cannot be spec- 
ified so that the observed cycloidal-plane flop cannot be 
explained. 

Below we summarize puzzling issues in the experimen- 
tal phase diagrams of i^MnOa: 



• An issue how the 6c-cycloidal order is stabilized 
in the strongly distorted materials despite the fact 
that in the perovskite manganites, the c axis is al- 
ways a hard magnetization axis. 

• A mechanism of the spin-cycloidal-plane (electric- 
polarization) flop with increasing GdFeOs-type dis- 
tortion. 

• A mechanism of the thermally induced cycloidal- 
plane flop or an issue why the 6c-cycloidal phase 
appears above the a6-cycloidal phase in tempera- 
ture. 

• Origin and nature of the sinusoidal coUinear phase 
in the intermediate temperature region. 

• An issue why the Mn spins in the AFM(A) and 
the sinusoidal collinear states direct along the or- 
thorhombic h axis [13, IM HI, IH, [H . 

These issues are interesting not only for potential appli- 
cations but also for fundamental physics. 

Previously, Ishihara and coworkers theoretically stud- 
ied the phase diagram with the use of a two-dimensional 
J1-J2 Ising model [sj. Although they predicted long 
wave-length magnetic structures between AFM(A) and 
AFM(E) phases, they inevitably failed to reproduce the 
cycloidal spin order and the sinusoidal collinear spin or- 
der since they treated the Mn spins as Ising spins. Re- 
cently, Dagotto and coworkers theoretically studied the 
phase diagram and the cycloidal spin order by employ- 
ing a two-orbital double-exchange model on the two- 
dimensional lattice with some additional terms pH . isH i . 
In Ref. they incorporated the lattice elastic energy 
as well as the DM interactions with vectors on the Mn- 
0-Mn bonds coupling to the oxygen displacements, but 
neglected the second-neighbor spin exchanges. They dis- 
cussed that the cycloidal spin order is stabilized due to 
the ferri-type arrangement of DM vectors realized by uni- 
form shifts of the oxygens. In Ref. [s^, they considered 
quite weak next-neighbor spin exchanges, but neglected 
the DM interactions. They reproduced the successive 
emergence of three magnetic phases in the ground state, 
i.e., AFM(A), cycloidal, and AFM(E) phases. However, 
in those studies, neither the spin-cycloidal-plane flop nor 
the sinusoidal collinear state was reproduced. In addi- 
tion, the issues on the easy-axis spin anisotropy in the 
AFM(A) phase and the emergence of 6c-cycloidal spin 
state were not addressed. Therefore, all of the above- 
listed issues still remain unresolved. 

In this paper, we theoretically investigate the mag- 
netoelectric phase diagrams of i?Mn03 in the absence 
of external magnetic held. We construct a microscopic 
spin model, which is basically a classical S'=2 Heisenberg 
model on a cubic lattice, but includes some additional 
magnetic anisotropics and interactions [36j . By analyz- 
ing this model using the Monte-Carlo method, we obtain 
phase diagrams in good agreement with experiments. We 
discuss a number of puzzling issues on the magnetoelec- 
tric phase diagrams of the perovskite manganites listed 
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above by particularly focusing on the mechanism of the 
spin-cycloidal-plane flop. We reveal that the DM inter- 
action between spins neighboring along the x or y axis 
as well as the single-ion anisotropies, because of which 
the c axis becomes a hard magnetization axis, favor the 
a6-cycloidal spin state with Pa, while the DM interaction 
between spins neighboring along the c axis favors the bc- 
cycloidal spin state with Pc. Their competition turns out 
to be controlled by the GdFeOs-type distortion, or equiv- 
alently by the second-neighbor exchanges J2 enhanced by 
this lattice distortion. This leads to a polarization flop 
from Pa to Pc with decreasing i?-site radius in agreement 
with the experiments. 

Here we note that although recent intensive experi- 
ments on the present manganites have uncovered rich 
phase diagrams and phenomena also in the presence of 
magnetic fields 



18, 33, M m, m, m, m m m, m 



\M\M,\M\M\M\^iMMM' these issues are beyond the 
scope of this paper. However, a microscopic model con- 
structed in this paper would offer a very powerful basis 
for studying the magnetic-field effects in the multiferroic 
manganites as well. 

Moreover, several fascinating phenomena originating 
from the strong magnetoelectric coupling have been ob- 
served in the multiferroic phases in /JMnOa, and have 
attracted appreciable interest as listed below: 

• The polarization flop can be induced also by an 
external magnetic field. In TbMnOa and DyMnOa, 
for instance, the electric polarization flops from Pc 
to Pa by app lying a magnetic field along the a or 
b axis [il, l33|. On the other hand, the polarization 
in Eui_a,Ya;Mn03 with x =0.4-0.5 flops from Pa to 
Pc with a magnetic field along the a axis [l^. In 
addition, control of the electric polarization vector 
using the rotating magnetic field was demonstrated 
experimentally ,38,., ,39] . 



• There are experimental [54 . l5a | and theoretical |56l . 
[s?! arguments on possible low-lying spin excita- 
tions activated by the electric-field component of 
the light (termed electromagnon) . Far- infrared 
spectroscopy observed corresponding light absorp- 
tions in the ferroelectric phases in TbMnOs ISa . 
m, M, DyMnOg GdMnOg [H, 
Gdi_:rTb:rMn03 ^ and Eui-^rY^MnOa [lllil 
at terahertz frequencies. 

• Giant magnetocapacitance was observed at a 
threshold ma gnet ic field for polarization flop in 
DyMnOg 0,113, [S, 

and its mechanism has 
attracted interest |4ll . [4^ . On the basis of 
the dielectric-dispersion measurements, this phe- 
nomenon was attributed to the electric-field-driven 
motion of the multiferroic domain wall between 
domains with orthogonal spin-cycloidal planes 
and concomitant orthogonal ferroelectric polariza- 
tions [4l|. It was revealed that the macroscopic 
motion of the multiferroic domain wall is possible 
with electric fields of practical magnitude, which 
enables the electric control of magnetic domains. 



It is widely recognized that interplay between charge, or- 
bital, spin and lattice degrees of freedom is a key to un- 
derstanding of these phenomena, similarly to the colossal 
magetoresistance in the hole-doped i?Mn03 [6^ . There- 
fore, construction of a microscopic model would represent 
a very important basis to approach these issues. Upon 
construction of a model, whether the experimental phase 
diagrams can be reproduced or not shall be an efficient 
test for its validity. 

Furthermore, new insights into the magnetoelectric 
coupling obtained here would be useful when we study 
the related phenomena not only for the present man- 
ganites but also for other newly-discovered multiferroic 
materials. The magnetically-induced ferroelectricity has 
recently been found in several transverse-spiral r nag - 
nets [63, lei such as NiaVaOg i69i], MnW04 [zO, iSf, 
LiCu202 III,!!!, LiCuV04 [Zl and CuO [zl^nd also 
in a transverse cone-spiral magnet GoCr204 f76l|. Fur- 
thermore, it was demonstrated that a conical spin order 
induced by an external magnetic field also generates the 
electric polarization in ZnGr2Se4 [T^I and Y-type hexa- 
ferrites [TSl . [79j . Under this circumstance, clarification 
of the essential physics behind the rich phase diagrams, 
and a microscopic description of the magnetoelectric sys- 
tems in the perovskite manganites have attained increas- 
ing importance because the manganite system is the first 
discovered example, and hence a prototype of a series of 
these new multiferroic materials. 

The organization of this paper is as follows. In Sec. [Ill 
we construct a microscopic spin model to describe the 
Mn 3c?-spin system in the perovskite manganites. In 
Sec. lIIIi we discuss the calculated results obtained by an- 
alyzing this model in the Monte-Carlo simulations. We 
study phase diagrams and properties of several phases 
in Sec. IIIIBI We address an issue on the spin-cycloidal- 
plane flop (or the electric-polarization flop) in Sec. IIII CI 
In Sec. IIII Dl we also investigate detailed ground-state 
spin structures of the cycloidal spin states by numeri- 
cally solving the Landau-Lifshitz-Gilbert equation. On 
this basis, we discuss the experimentally claimed ellipti- 
cal modulation of the spin cycloid [10., Ml]- In Sec. IIII El 
by examining effects of the single-ion anisotropies, we ad- 
dress an issue why spins in the AFM(A) state and those 
in the sinusoidal coUinear state are aligned parallel to 
the b axis. In Sec. IIIIFI we discuss the nature and ori- 
gin of the sinusoidal coUinear spin phase at intermediate 
temperatures. Section IIVI is devoted to conclusion and 
discussion. We also explain how we calculate values of 
parameters in our model in Appendices A, B and C. 



II. MODEL AND METHOD 

A. Model Hamiltonian 

To describe the Mn 3(i-spin system in /JMnOa, we em- 
ploy a classical Heisenberg model with some additional 
interactions and magnetic anisotropies on a cubic lattice. 
In this Hamiltonian, we treat the Mn 5=2 spins as clas- 



sical vectors, i.e. Si={yJ S"^ — S^cosOi, yj S'^ — S'^ sin^i, 
Sc) with respect to the Pbnm o,, b and c axes. The Hamil- 
tonian consists of four contributions as 



with 



Tic 



(2) 



(3) 



TiuM 

Ticuh 



— Jab Si ■ Sj + J2 Si ■ Sj 

<i,j> <i,j> 
c 

+Jc ^ ^ Si ■ Sj, 

<hj> 

Dj2sl + Ej2i-iy^+^y{Sl^Sl), (4) 

i i 

d% ■ {S, X S,), (5) 

<i,j> 

^^^^TjTYy + 'S'yi + 'S'^i), (6) 



where ix, iy and iz represent coordinates of i-th Mn ion 
with respect to the cubic x, y and z axes. In the foUow- 
ing, we explain each of these terms in detail, including 
definitions of other indices and variables. 



B. Superexchange interactions 

The first term Ticx represents superexchange interac- 
tions. This term contains ferromagnetic exchanges Jat 
on the Mn-Mn bonds along the x and y axes, antiferro- 
magnetic exchanges J2 on the in-plane diagonal Mn-Mn 
bonds along the orthorhombic b axis, and antiferromag- 
netic exchanges Jc on the Mn-Mn bonds along the c axis 
see Fig.inja). 

The Jahn- Teller distortion in i?Mn03 causes a C-type 
orbital ordering. Consequently, the ferromagnetic Jab is 
caused by the antiferro arrangement of eg orbitals in the 
ab plane, while along the c axis the antiferromagnetic 
Jc is caused by the ferro-orbital stacking. On the other 
hand, the antiferromagnetic J2 arises from an exchange 
path between two Mn eg orbitals in the b direction via two 
oxygen 2p orbitals as shown by a dashed hue in Fig.[3fb). 
The strength of the second-neighbor exchanges J2 in- 
creases with increasing GdFeOs-type distortion since the 
distortion enhances hybridization between these two oxy- 
gen 2p orbitals [1 [sif ■ 

It should be mentioned that in other orthorhombic per- 
ovskite compounds, the strength of the second-neighbor 
exchanges J2 is negligibly weak as compared to that of 
the nearest-neighbor exchanges Jab- By contrast, in the 
manganites, these two exchanges J2 and Jab are compa- 
rable in strength, and the second-neighbor exchanges J2 
play important roles in determining the magnetic prop- 
erties. This can be attributed to the electronic struc- 
ture of Mn'^+ ion with ^29^9 electron configuration. For 
the nearest-neighbor Mn-Mn pairs, there are two oppo- 
site contributions to the superexchange Jab as shown in 




(b) 




FIG. 3: (Color online) (a) Superexchange interactions in 
RMnOs and tilted local coordinate axes ^i, rji, and ("i attached 
to the i-th MnOg octahedron. With regard to the superex- 
change interactions, we consider ferromagnetic exchanges Jab 
on the Mn-Mn bonds along the cubic x and y axes, anti- 
ferromagnetic exchanges J2 on the in-plane diagonal Mn-Mn 
bonds along the orthorhombic b axis, and antiferromagnetic 
exchanges Jc on the Mn-Mn bonds along the c axis, (b) Con- 
figuration of the occupied Cg orbitals and superexchange in- 
teractions {Jab and Jc, solid curves) in the ab plane. Exchange 
path for the superexchange J2 via two oxygen 2p orbitals is 
indicated by a dashed line. 



(a) 



3y2^ 



t2g 



(b) 



3x2-r2 



t2g 



3y2-r2' 



t2g 



Jab = J^M - j'^^Af 



J2=J^»AF + /*AF 



3y2-r- 



t2g 



FIG. 4: (Color online) (a) Opposite contributions to the 
nearest-neighbor ferromagnetic exchange Jab from the t2g- 
and Cg-orbital sectors, (b) Cooperative contributions to the 
second-neighbor antiferromagnetic exchange J2 from the t2g- 
and Cg-orbital sectors. 
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TABLE I: Calculated superexchange parameters [Jat and Jc) 
and single-ion anisotropy parameters {D and E) for several 
iiMnOs compounds. Structural data used in the calculations 
are taken from Refs. [84l . Issl . |86| , [stI . IssI ]. In the second row, 
their ground-state magnetic structures are presented. Cal- 
culated values of D and E contain some ambiguities due to 
uncertainty in the optical absorption data [89!, ydl . The val- 
ues in the upper (lower) rows are calculated using the data 
from Ref. }89j (Ref. Values of the second-neighbor anti- 

ferromagnetic exchange J2 are also presented in the last row, 
which are calculated using the J1-J2 model, (see Apendix B) 
The value of J2 in EuMnOs [AFM(A)] is deduced by extrap- 
olation. 



RMnOa 


EuMnOa 


TbMnOa 


DyMnOs 


YMnOs 


Magnetic state 


AFM(A) 6c-cycloidal fec-cycloidal 


cycloidal? 


Jab (meV) 


0.85 


0.79 


0.76 


0.77 


Jc (meV) 


1.27 


1.26 


1.20 


1.25 


D (meV) 


0.32 


0.33 


0.33 


0.34 




0.24 


0.25 


0.25 


0.25 


E (meV) 


0.34 


0.34 


0.34 


0.34 




0.25 


0.24 


0.24 


0.24 


J2 (meV) 


~0.5 


0.62 


0.92 


1.9 



Fig. mja). The exchange between the 5=1/2 spins on 
the staggered eg orbitals gives a ferromagnetic contribu- 
tion, while in the t2g-orbital sector, the antiferromagnetic 
couphng is reahzed between the S—3/2 spins. The can- 
cellation of these opposite contributions results in the 
small magnitude of Jab- On the other hand, for the 
Mn-Mn pairs neighboring along the b axis, both the Cg- 
orbital and the t2g-orbital sectors give antiferromagnetic 
contributions to the superexchange J2 as in Fig. HJ^b), 
resulting in the relatively large magnitude of J2. The 
frustration between the ferromagnetic Jab and the anti- 
ferromagnetic J2 stabilizes a spiral spin order in i?Mn03 
compounds with relatively strong GdFeOa-type distor- 
tion like TbMnOa, DyMnOs, and certain kinds of solid 
solutions. 



We calculate the values of superexchange parameters 
Jab and Jc for several i?Mn03 compounds by using for- 
mulae given in Refs. [H, [s^. The results are listed in 
Table HI The structural data used in this calculation are 
taken from Refs. M, 85., M, iSS, [11. The changes of 
their values are very small upon the i?-site variation as 
far as vicinities of the multiferroic phases are concerned. 
We take Jab=0.80 meV and Jc— 1.25 meV in the follow- 
ing calculations. We also estimate the values of J2 in 
TbMnOa, DyMnOa, YMnOa and several sohd solutions 
by employing the two-dimensional J1-J2 classical Heisen- 
berg model, which reveals that the value of J2 monoton- 
ically increases with decreasing i?-sitc radius. For the 
details, see Appendices A and B. 



TABLE IL Structural parameters of EuMnOa from Ref. [871. 
a (A) b (A) c (A) Xq^ yo^ Xq^ yp^ zp^ 
5.3437 5.8361 7.46186 0.0974 0.4714 0.7055 0.3247 0.04845 



C. Single-ion anisotropy 

The second term Tisia denotes the single-ion 
anisotropics, which is determined by the wave function of 
occupied eg orbital, or equivalently by the local environ- 
ment of Mn^+ ion surrounded by six oxygens. Here ^i, 
1)1 and Ci are tilted local axes attached to the i-th MnOg 
octahedron as shown Fig.[3Ua). This term consists of two 
parts, i.e. H^g^ and H^g^. The former part implies 
that the C,i axis is a local hard magnetization axis at every 
site, and consequently the crystal c axis becomes a hard 
magnetization axis since the Q axis at each site directs 
nearly along the c. On the other hand, the latter part 
■^^g implies that the and m axes become a local hard 
magnetization axis alternately in the ah plane. This is 
due to the staggered ordering of occupied eg orbitals, i.e. 
the 3a;^-r^/3j/^-r^ type orbital ordering. 

We derive directional vectors of the , r\i and axes 
with respect to the a, h and c axes using coordination 
parameters of oxygens as 



a[0.25 + (0.75 - J] 

6[0.25-(-l)'-+*«(yo. -0.25)] 

a[-0.25 + (-l)^-+*« (0.75 - xo^) 
6[0.25+(-l)'-+'='(yo. -0.25)] 



-c(-l)*-+^"+'-zo2 



C = 



— a(— 1)*="+*"+*^ xoi 
&(-l)'^(0.5-yoJ 
0.25c 



(7) 



(8) 



(9) 



Here , 2/O2 and are the coordination parameters 
of the in-plane oxygens, a;oi , J/Oi and zq^ are those of the 
out-of-plane oxygens, and a, h and c are the lattice pa- 
rameters. For values of these parameters, we use the ex- 
perimental data of EuMnOa 8^ (see Table [Til) through- 
out the calculations for simplicity. We have confirmed 
that the results are not significantly changed even if we 
use the values for other iJMnOs compounds especially 
with i?=Gd, Tb, Dy, Y and Ho. We note that the axis 
vectors r]i and C,i in Eqs. Q-® are not normahzed. 
Using these vectors, we can rewrite H.^^^ and as 

n^^ = i?5](5,.c/K.|)', (10) 

i 

= i?5](-i)^==+^« (11) 

i 

X \{s,■i^l%\f-{s,■^.|\r].\f\. 

In Table. [H the values of D and E are listed, which are 
calculated using formulae given in Ref. [11]. Noticeably, 
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djk = (-ac, -Pc 0) 



dsa = (-ac, Pc 0) 




The five DM parameters, Uab, Pat, lab, cic, and /3c, 
in LaMnOa was evaluated in the first-principles calcula- 
tion (9JI , which shows that the DM vector on the in-plane 



I Q dn = {-aab, Pat, yab) 
d'A = (aafe, Pafe, yab) 



FIG. 5: (Color online) Dzyaloshinsky-Moriya vectors associ- 
ated with different Mn-O-Mn bonds, which are expressed by 
five parameters, aab, Pat, 'yab, olc and /3c- The vectors on the 
in-plane Mn-O-Mn bonds direct nearly perpendicular to the 
plane consisting of the corresponding Mn, O and Mn ions. 



the values are all approximately equal in the compounds 
located near/in the multifcrroic phases in the phase dia- 
gram although the calculated values contain some am- 
biguities due to uncertainty in the optical absorption 
data pl,l90l|. We take D^^.lh meV and -B^O.SO meV in 
the following calculations. For details of the parameter 
calculations here, see Appendix C. 



D. Dzyaloshinsky-Moriya interaction 

The third term TYdm denotes the Dzyaloshinsky- 
Moriya (DM) interactions [9l|, [sllll. The DM vector 
rfg is defined on the Mn(i)-0-Mn(j) bond along the a 
axis (a=x, y and z), and the following antisymmetric re- 
lation holds; — —df^ Given that the Mn(j) ion is 
adjacent to the Mn(i) ion in the positive a direction, the 
DM vectors are given by (see also Fig. [5]) 



Mn-O-Mn bond is nearly perpendicular to the plane con- 
sisting of the corresponding Mn, O, and Mn ions, while 
the DM vector on the out-of-plane bond is not. This can 
be understood as follows. The spin exchange on the in- 
plane Mn-Mn bond is governed by the exchange path via 
the inbetween one oxygen, i.e. the Mn-O-Mn path. Thus 
the direction of the DM vector on the in-plane bond is 
dominantly determined by local symmetry of the Mn-O- 
Mn bond. On the other hand, the spin exchange on the 
out-of-plane Mn-Mn bond is not necessarily governed by 
the Mn-O-Mn path along the c axis, and contributions 
from indirect paths containing more than one oxygens are 
not negligible since the occupied eg orbitals are directed 
not along the c axis but in the ab plane. As a result, the 
DM vector on the out-of-plane bond does not reflect the 
local symmetry of the Mn-O-Mn bond. 

In addition, according to the electron-spin resonance 
(ESR) measurements for LaMnOs [oi, HUl, the DM 
vector on the out-of-plane bond is approximately four 
times larger in magnitude than the vector on the in-plane 
bond. We expect that the above hndings also hold in 
other iJMnOs compounds although neither theoretical 
nor experimental studies on the DM vectors presently 
exists except for i?=La. In our calculation, we take 
aafc=0.10 meV, /5afc=0.10 meV, 7ab=0.14 meV, ^^=0.30 
meV and /3c=0.30 meV. With this set of parameters, the 
above two features are reproduced: the DM vector on the 
in-plane bond is nearly perpendicular to the correspond- 
ing Mn-O-Mn plane, and the vector on the out-of-plane 
bond is approximately three times larger than the vector 
on the in-plane bond. 

We also mention that the a components of DM vectors 
on the out-of-plane bonds have the same signs within 
each plane, but the signs alternate along the c axis. This 
situation gives rise to a weak ferromagnetism with mo- 
ments along the c axis due to the spin canting in the 
AFM(A) phase in agreement with the experimental ob- 
servation [H, [oil . 



dl 



d^ 



(-l)'^'+*«7afc 

-{-ly^+'y+'-aab 

i-lY^+'-lab 





(12) 



(13) 



(14) 



because of the crystal symmetry. 



E. Cubic anisotropy 

The last term Ticub represents the cubic anisotropy, 
which comes from nearly cubic symmetry of the per- 
ovskite lattice. Here, Sx, Sy and Sz are written as 
Sx = l/V2iSb - Sa), Sy = l/V2iSb + Sa), and = Sc. 
The coupling constant a in the Mn3+ ion surrounded by 
the octahedrally coordinated oxygens was evaluated to 
be 0.0162 meV in the ESR measurement [131 . We ne- 
glect a slight contribution from the orthorhombic lattice 
distortion to the spin anisotropy since it is expected to 
be very small. 
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TABLE III: Model parameters used in the calculations for 
each term of the Hamiltonian. The energy unit is meV. 



Hex 


Jab = O.SO, Jc-- 


=1.25 






75=0.25, 


=0.30 




Ti-oM 


aab=0.10, Pab 


=0.10, -yab 


=0.14 






=0.30, 7c 


=0.0 


Ti-cub 


0=0.0162 







F. Method 

We calculate thermodynamic properties of the model 
([2|) by using the Monte-Carlo method. The model pa- 
rameters used in the calculations are summarized in Ta- 
ble [Till To avoid a critical slowing down in the frustrated 
systems, we employ the replica exchange Monte-Carlo 
method "OO*!. We take an exchange sampling after ev- 
ery 400 standard Monte-Carlo steps. Typically, we per- 
form 600 exchanges after sufficient thermalization Monte- 
Carlo steps. In the following, we mainly show the results 
obtained for systems with 48 x 48 x 6 sites under the peri- 
odic boundary condition. By performing the calculations 
also for systems with 36x36x6 sites at some parameter 
values, and by adopting the open boundary condition 
also, we confirm that the finite-size effect is small enough 
and never affects our conclusion. 

On the other hand, we study ground-state properties 
by numerically solving the Landau-Lifshitz-Gilbert equa- 
tion. We derive effective local magnetic fields H^^ acting 
on the i-th Mn spin Si from the Hamiltonian Ti as 



Hf = -dH/dS,. 



(15) 



Then, we construct a Landau-Lifshitz-Gilbert equation 
with thus obtained local fields; 



dt 



^s, X h: 



cff 



s 



where ao is the dimensionless Gilbert-damping coeffi- 
cient introduced phenomenologically. For the value of 
ao, we take a rather small value of ac — 0.01 to achieve 
a slow relaxation towards a stable spin structure with 
a minimum energy. We solve this equation using the 
Runge-Kutta method. For the convergence, we use ther- 
mally relaxed spin configurations obtained in the Monte- 
Carlo simulations at the lowest temperature (/cbT ~ 0.5 
meV) as initial states. 



III. RESULTS 
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FIG. 6: (Color online) Theoretically obtained magnetoelec- 
tric phase diagram in plane of temperature and J2. Here, 
AFM(A)+WFM denotes the A-type antiferromagnetic phase 
with weak ferromagnetism due to the spin canting, and PE 
and FE denote paraelectric and ferroelectric phases expected 
in the spin-current model [9I [Til. [T^. respectively. 



Under this assumption, we investigate the T-rn diagrams 
for i?Mn03 (r^j is the ionic i?-site radius) and the T- 
X diagrams for solid solutions by varying the value of 
J2. In Sec IIIIBl we first display a theoretically ob- 
tained T-J2 phase diagram, which shows good agree- 
ment with the experimentally obtained T-x phase dia- 
gram of Eui_a;Ya;Mn03 in Ref. [l^. Then we show cal- 
culated results for the momentum dependence of spin 
and spin-helicity correlation functions, and the tempera- 
ture dependence of specific heat and spin-helicity vector. 
We discuss how we identify phase transitions and mag- 
netic structures in the phase diagram with the aid of 
these results. In Sec. IllICi we discuss a mechanism of 
the orthorhombic-distortion-induced spin-cycloidal-plane 
(electric-polarization) flop. We also demonstrate that the 
regime of 6c-cycloidal spin (P||c) phase increases as the 
DM parameter ac increases, and the experimental phase 
diagram of the Gdi_a;Tba;Mn03 system in Ref. [l^ is re- 
produced well for a rather large value of ac= 0.38 meV. 
In Sec. IIIIDl we study the ground-state spin structures 
in the ab- and 6c-cycloidal spin phases, and discuss the 
experimentally claimed elliptical modulation of the cy- 
cloidal spin structure. In Sec. IIII El by examining effects 
of the single-ion anisotropies, we address the issue why 
the Mn spins are aligned along the b axis in the AFM(A) 
and sinusoidal collinear phases. In Sec. IIIIFl we discuss 
the sinusoidal collinear spin phase in the intermediate 
temperature region by examining the role of the single- 
ion anisotropies on its emergence. 



A. Overview 

When we study the experimental phase diagrams in- 
troduced in Sec. [H we assume that main roles of the 
GdFeOs-type distortion on the magnetic and electric 
properties in i?Mn03 are induction and enhancement 
of the second-neighbor antiferromagnetic exchanges J2. 



B. Phase diagram 

In Fig. [6l we display a theoretically obtained T-J2 
phase diagram, which reproduces the experimental T-x 
phase diagram of Eui_2:Y2;Mn03 system [Fig.[5][b)]. This 
theoretical diagram shows that at the lowest tempera- 
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ture, the canted AFM(A) state emerges for J2 < 0.59, 
while for J2 > 0.59, the cycfoidal spin states emerge. 
For 0.59< J2 < 0.64, the system undergoes two ther- 
mal phase transitions, and the paramagnetic, sinusoidal 
coUinear, and a6-cycloidal phases successively emerge 
with lowering temperature. On the other hand, for 
J2 > 0.64, the system undergoes three thermal phase 
transitions, and four magnetic phases, i.e., the paramag- 
netic, sinusoidal coUinear, a6-cycloidal and 6c-cycloidal 
phases, successively emerge as temperature decreases. 
The third transition is nothing but a spin-cycloidal-plane 
flop from the be plane to the ab plane. The regime of bc- 
cycloidal phase increases with increasing J2. According 
to the spin-current model, the ferroelectricity with Pa is 
expected to show up in the a6-cycloidal spin phase and 
that with Pc in the &c-cycloidal spin phase. Note that 
in the gray region in Fig. [6l we obtain several compli- 
cated magnetic structures, whose spin correlation func- 
tions show peaks at several k-points. They are considered 
to be artifacts of the finite-size calculation for 48x48x6 
sites. In the vicinity of the AFM(A) phase, pitches of the 
spiral and sinusoidal spin states become so long and the 
magnetic unit cells become so large that the calculation 
tends to be seriously affected by the finite-size effects. 

We identify each of the phases and each of the transi- 
tions in the above phase diagram by calculating several 
physical quantities. The magnetic structures are assigned 
from calculated momentum dependence of the spin cor- 
relation function Sa{k,T) and that of the spin-helicity 
correlation function H^{k,T). On the other hand, the 
transition points are determined from calculated temper- 
ature dependence of the specific heat Cs{T) and that of 
the spin-helicity vector hll{T). 

The spin and spin-helicity correlation functions 
Sa{k,T) and H^{k,T) for a=a, b, c are calculated by 

S„{k,T) = (17) 



Hi{k,T) = _^(/,^,/,^^.)e*(— (18) 

where the bracket denotes the thermal average. Here, 
is the a component of local spin-helicity vector h\ — 
i^ai^ ^bi^ ^ci): which is defined as 

h'^ = {S, X S,+t)/S^. (19) 

In the following, we write these correlation functions sim- 
ply as Sa{k) and H'^{k) by omitting T. 

In a cycloidal spin state, all of the local spin-helicity 
vectors point in the same direction perpendicular to the 
basal cycloidal plane. In this sense, the cycloidal spin 
state can be regarded as a ferrohelicity state. Further- 
more, the a6-cycloidal spin state has a ferrohelicity com- 
ponent along the c axis, while the &c-cycloidal spin state 
has that along the a axis. Thus the spin-helicity c- 
component correlation H^{k) has a peak at fc=(0, 0, 0) in 



the a6-cycloidal spin phase, while in the 6c-cycloidal spin 
phase, its a-component correlation H^{k) has a peak also 
at fc=(0, 0, 0). In the following, we overview features of 
Sa{k) and Ha{k) in each magnetic phase by taking typ- 
ical points in the phase diagram as examples, which are 
indicated by cross symbols in the inset of Fig.[7l^a). 

In Fig. [TJa) we display the calculated Sa{kx, ky,TT) 
(upper left panel) and H^{kx,ky,0) (right panel) in the 
AFM(A)-|-WFM phase. In this phase, the Mn spins are 
aligned nearly parallel to the b axis, and couple ferromag- 
netically in the ab plane with antiferromagnetic stacking 
along the c axis. As a result, the spin 6-component cor- 
relation Sb{k) has a sharp peak at fe = (0,0, tt), while 
Sa{k) and Sc{k) do not. On the other hand, the spin- 
helicity correlation Ha{k) has no structure for any a. In 
addition, weak ferromagnetic moments arise due to the 
spin canting. As shown in the lower panel of Fig. [Tj^a), 
the spin c-component correlation Sc(k) exhibits a tiny 
but sharp peak at fc = 0. Moreover, this spin canting in 
fact generates finite local spin helicitics h': between spin 
pairs neighboring along the c axis, defined as 

/i^ = {S, X S,+c)/S\ (20) 

These vectors direct in the positive a direction or in the 
negative a direction uniformly in the plane, and these 
two kinds of planes are stacked alternately. Thus if we 
introduce the correlation function of h'^ as 

H^^{k,T) = -l^(/i^^/i=^^.)e*("--.), (21) 

the a-component correlation H^{k,T) has a peak at 
fc=(0, 0, tt) not shown. 

In Fig. [ZUb), we display the calculated Sa{kx, ky,7T) 
and H^{kx, ky, 0) in the afo-cycloidal phase. In this phase, 
the Mn spins rotate in the ab plane, while they couple 
antiferromagnetically along the c axis. As a result, the 
spin correlations Sa{k) and S'h(fe) exhibit sharp peaks at 
{ka,kb,kc) — (0,±(3'm,7r) with a finite q„i. The rotating 
spins in the ab plane give rise to a ferro-arrangement of 
the spin helicities along the c axis. This results in a 
peak of spin-helicity correlation H^ik) at A;=0. 

The correlation functions in the &c-cycloidal phase in 
Fig- Eljc) exhibit that the spin correlations Sh{k) and 
Sc{k) show sharp peaks at (fc^, kb^ kc) = (0, ±gm, tt), and 
the spin-helicity correlation H^{k) has a sharp peak at 
fc = 0. This is because the Mn spins rotating in the be 
plane give rise to a ferro-arrangement of the spin helici- 
ties hll along the a axis. 

Finally, the results for sinusoidal coUinear phase are 
shown in Fig. [7^d). In this phase, the coUinearly aligned 
Mn spins (||&) are sinusoidally modulated with a finite 
modulation vector along the b axis. This results in 
sharp peaks of the spin 6-component correlation Sb{k) 
at {ka,ki,,kc) = (0, ±qm,7r) with a finite q™, as well as 
no remarkable structure in the spin-helicity correlations 
H^{k) for any a. 
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(a) AFM(A)+WFM phase (J2=0.44 meV, /(b7=0.5 meV) 
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(b) sinusoidal collinear spin phase {J2=0.80 meV, ksT =3.5 meV) 
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(c) 5c-cycloidal spin phase {J2=0.80 meV, ksT =2.5 meV) 
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(d) at)-cycloidal spin phase (J2=0.80 meV, ksT =0.5 meV) 
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FIG. 7: (Color online) Calculated spin correlation functions Sa{k) for kc='K (left panels) and spin-helicity correlation functions 

Ha(k) for fcc=0 (right panels) in the k^-ky plane {a=a, b and c) for different magnetic states at the points indicated by 

cross symbols in the inset; (a) the canted AFM(A) state [AFM(A)+WFM state] at J2=0.44 meV and fcBT=0.5 meV, (b) the 
sinusoidal collinear spin state at J2=0.80 meV and fcBr=3.5 meV, (c) the fec-cycloidal spin state at J2=0.80 meV and fcBT=2.5 
meV, and (d) the a6-cycloidal spin state at J2=0.80 meV and ^3^=0. 5 meV. For the AFM(A)+WFM state, spin correlation 
function Sa{k) for kc—O is also displayed, which indicates small ferromagnetic moments along the c axis. 
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FIG. 8: (Color online) Calculated temperature profiles of spe- 
cific heat Cs (T) and spin-helicity vector (T) for various val- 
ues of J2 (along the vertical dashed lines in inset); (a) J2=0.36 
meV, (b) J2=0.62 meV, and (c) J2=0.80 meV. Here, h.y{T) 
(7=0, b, c) denotes the 7 component of h'^{T). In the ab- 
cycloidal [6c-cycloidal] phase, the value of hc{T) [ha{T)] is 
much larger than other two components. 



In Fig. [51 we show the temperature dependence of spe- 
cific heat Cs (T) and spin-helicity vector /i^ (T) for various 
values of J2, which are calculated by, 



Cs{T) 



d{n)/d{k^T), (22) 
U\Y^S,^S.+,\)IS\ (23) 



N 
1 

iV 



Here, hj(T) {j—a, b and c) denotes the 7 component 
of spin-helicity vector h.^(r) = [ha{T), hb(T), hdT)]. In 
the a6-cycloidal [&c-cycloidal] phase, its c [a] component 
hc{T) [ha{T)] has a large value, while other two compo- 
nents are strongly suppressed. On the other hand, in the 
paramagnetic, canted AFM(A) and sinusoidal coUinear 
phases, all of these three components should be nearly 



equal to zero. 

In Fig. [Hta), we can see a single phase transition in 
Cs{T) for a small value of J2=0.36 meV, which corre- 
sponds to a transition from paramagnetic into canted 
AFM(A) phases. The spin-helicity vector hll{T) is nearly 
equal to zero constantly through this transition. 

For a larger value of J2=0.62 meV, Cs(T) shows two 
peaks as shown in Fig. [Sfb), indicating that successive 
two transitions occur with lowering temperature. The 
spin-helicity vector hll{T) is almost unchanged, and is 
approximately zero through the first transition, at which 
the system enters into the sinusoidal collinear phase from 
paramagnetic phase. Contrastingly, at the subsequent 
transition, its c component hc(T) starts increasing, indi- 
cating a transition into the afe-cycloidal phase. 

For a further increased value of J2=0.80 meV, the sys- 
tem exhibits successive three phase transitions as temper- 
ature decreases see Fig.[5{c). The first one is a tran- 
sition from paramagnetic to sinusoidal collinear phases 
through which all of the three components of h^{T) are 
again approximately zero constantly. At the second tran- 
sition, its a component, ha{T), increases, while other two 
components are remain to be small, indicating a transi- 
tion into the &c-cycloidal phase. With further lowering 
temperature, the a component ha(T) suddenly drops, 
while the c component hc{T) steeply increases at the 
third transition point, indicating a spin-cycloidal-plane 
flop from the be plane to the ab plane. 



C. Orthorhombic-distortion-induced polarization 
flop 

In iJMnOs, the c axis is always a hard axis for magneti- 
zation due to the single-ion anisotropics Ti^^- Therefore, 
the &c-cycloidal spin state seems to be higher in energy 
than the a6-cycloidal spin state at first sight. In this 
sense, the emergence of &c-cycloidal spin order is puz- 
zling. Moreover, the 6c-cycloidal spin phase emerges next 
to the a6-cycloidal spin phase with increasing GdFeOa- 
type distortion. The mechanism of this lattice-distortion- 
induced cycloidal-plane flop has not been clarified yet. 

If the Hamiltonian contains the superexchange term 
only, spins in the cycloidal state rotate uniformly. In the 
ground state, the angle 4> between adjacent two spins 
along the x ot y axis satisfies the relation, cos (f> = 
Jab/(2J2). On the other hand, if we introduce the DM 
interactions, the rotation angles should be modulated, 
and are no longer the same. 

In the a6-cycloidal spin state, the rotating spins cou- 
ple dominantly to the c components of DM vectors on 
the in-plane Mn-O-Mn bonds. Their magnitudes are all 
equal to jab, and their signs (i.e. +jab and —jab) alter- 
nate along the x and y bonds see the left panel of 

Fig-UJa). Without DM interaction, the spins rotate with 
the same rotation angles of (pab- On the other hand, in 
the presence of DM interactions, the rotation angles be- 
come to be alternately modulated into (f>ab + '^(f'ab and 
(pab — ^<f>ab with A<j)ab > to get an energy gain from the 
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FIG. 9: (Color online) (a) [(b)] Left panel: Spin structure in 
the a?)-cycloidal [bc-cycloidal] state and arrangement of the c 
[a] components of DM vectors on the in-plane [out-of-plane] 
Mn-O-Mn bonds. The symbols and express their signs, 
i.e. for the positive sign and for the negative sign. Right 
panel: (upper panel) Spin directions in the afc-cycloidal [bc- 
cycloidal] state with uniform rotation angles in the absence 
of DM interactions, and (lower panel) spin directions with 
modulated rotation angles in the presence of DM interactions. 

DM interactions see the right panel of Fig. O^a). We 

can derive the energy gain due to this angle modulation 
as 

AE^ll/N = -7afc52|sin(<?!)ah- A(/)a6)-sin(/)afc| 

= -TafcS'^ I cos 0ab I A(/)afc. (24) 

This expression implies that the energy gain jAE'g^l 
is reduced with increasing (pat because the prefactor 
|cos0ab| becomes maximum (=1) for (pat =0 but de- 
creases as increases. Thus the a6-cycloidal spin 
state is destabilized with increasing J2 or with increasing 
GdFeOs-type distortion. 

On the other hand, the spins in the &c-cycloidal state 
dominantly couple to the a components of DM vectors 
on the out-of-plane Mn-O-Mn bonds. Their magnitudes 
are all equal to ac, and their signs are the same within a 

plane, but alternate along the c axis see the left panel 

of Fig.[ni^b). Without DM interaction, the angles between 
adjacent two spins along the c axis are uniformly (pc — 
because of the strong antiferromagnetic coupling Jc- In 
the presence of DM interactions, the angles again suffer 
from staggered modulation into vr+A^c and tt— A0c with 

A0C > see the right panel of Fig.[5{b). Similarly to 

the a6-cycloidal case, we can derive the energy gain due 



to the angle modulation as 

I^E^m/N = -ac52|cos0,|A0, 

= -a,S''A(Pc. (25) 

Since the value of |coS(/)c| becomes maximum (=1) at 
(pc— TT, the energy gain AE^y^ for the &c-cycloidal spin 
state is always maximum irrespective of the value of J2 . 

As a result, the energetical advantage of 6c-cycloidal 
spin state relative to the a5-cycloidal one due to the DM 
interactions, |Ai?Q'jy[ — AiJ^^^j, increases as J2 increases. 
The 6c-cycloidal spin state is expected to be stabilized 
when the above energy difference dominates over the en- 
ergetical disadvantage due to the hard magnetization c 
axis. 

According to the above discussion, we expect that the 
a component of DM vector on the out-of-plane Mn-O-Mn 
bond is relevant to the stability of 6c-cycloidal spin state. 
To confirm this, we investigate T-J2 phase diagrams for 
various values of ac- Figures fTUTal- fTUf c) . together with 
Fig. [21 show that the regime of &c-cycloidal phase in- 
creases as ac increases. Note that the phase diagram for 
ac=0.30 meV depicted in Fig. [5] should be placed be- 
tween Figs. [TUT a) and fTUf b). which are for ac=0.24 meV 
and ac=0.34 meV, respectively. 

The experimental phase diagram of Gdi_a;Tba;Mn03 
depicted in Fig.[2)Jb) is reproduced when Q!c=0.38 meV as 
in Fig. [TCTc). This value is slightly larger than ac=0.30 
meV, for which the phase diagram of Eui_:rY2;Mn03 
is reproduced. This difference may be due to different 
i?-site-radius dependence of the GdFeOs-type distortion. 
According to Ref . , the lattice parameters of TbMnOa 
are equivalent to those of Eui_a;Y2:Mn03 with x ^ 0.85 
although TbMnOs is expected to be located at x ~ 0.4 
in the T-x diagram in terms of the averaged i?-site ra- 
dius. This indicates that the lattice of TbMnOa is more 
significantly distorted than expected from comparison to 

the Eui_a;Ya;Mn03 system we may have to consider 

not o nly t he average of i?-site radii but also their vari- 
ance |lO0j |. The out-of-plane Mn-O-Mn bond angles in 
Gdi_a;Tba;Mn03 tend to more significantly deviate from 
180° than those in Eui_a:Ya:Mn03 system, resulting in 
larger DM vectors on the out-of-plane bonds or a larger 
value of ac- (Note that the DM vector becomes zero 
when the bond angle is 180° because of symmetry.) 

We should also note that there is a slight difference 
between the experimental diagram of Gdi_a;Tba;Mn03 

and theoretical diagram for ac=0.38 meV compare 

Fig. [DJb) and Fig. [TUf c). In the experimental one, the 
phase boundary between the ah- and 5c-cycloidal spin 
phases slightly bends, and in the very narrow region 
near the phase boundary, the system exhibits a reen- 
trant behavior with successive transitions from he- to 
ah- and again to 6c-cycloidal spin phases with lowering 
temperature. In addition, the afo-cycloidal spin order is 
absent in the ground state. These points are not repro- 
duced in our calculation. This discrepancy may be solved 
by considering effects of the /-electron moments on the 
rare-earth ions, which order approximately below 10 K. 
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FIG. 10: (Color online) Theoretically obtained magnetoelec- 
tric phase diagrams for (a) ac=0.24 meV, (b) ac=0.34 meV, 
and (c) ac=0.38 meV. Here expresses magnitude of the a 
component of DM vector on the out-of-plane Mn-O-Mn bond. 
PE and FE denote paraelectric and ferroelectric phases ex- 
pected in the spin-current model [ol. [Til, [l^ . respectively. 



Here we emphasize that overall features of the phase di- 
agram of Gdi_2;Tba;Mn03 can be reproduced within the 
Mn 3(i-spin sublattice model. However, to reproduce its 
subtle features, we may need to consider couplings be- 
tween the Mn 3d spins and t he rare-earth /-electron mo- 
ments [i,[l0lliM[iQl,[l03- By contrast, in the case of 
Eui-ajYj^MnOs without interference from /-electron mo- 
ments, the agreement between the experiment and the 
calculation is quite good. 

To summarize this section, the ab- and 6c-cycloidal 
spin states in RMnOs are stabilized by the single-ion 
anisotropy or the DM interaction. Namely, the hard mag- 
netization c axis due to the single-ion anisotropy gives a 
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FIG. 11: Temperature dependence of y Sc{q)/Sb{q) in 
TbMnOs measured in the spin-polarized neutron-scattering 
experiment from Ref. [l^. Here q is the propagation wave 
vector of the 6c-cycloidal spin state. 



relative stability to the a6-cycloidal state, while the bc- 
cycloidal state is stabilized by the DM vectors on the out- 
of-plane bonds. On the other hand, the ac-plane spin cy- 
cloid is unfavorable. Recent first-principles calculations 
for TbMnOa also confirmed this tendency [H, [l3| • 



D. Magnetic structures in the cycloidal spin states 

The calculated spin correlation functions Sa{k) and 
Sb{k) in the a6-cycloidal spin phase are not equal in mag- 
nitude of the peak see Fig. ^d) . This is also the case 

for 5f,(fc) and Sc{k) in the fec-cycloidal spin phase see 

Fig.[5]^c). At first sight, this seems to mean that the mag- 
nitude of ordered magnetic moment depends on its direc- 
tion. The inequivalency between Sb{q) and Sc{q) was ac- 
tually observed in the recent spin-polarized neutron scat- 
tering experiments for TbMnOa and Tbi-^jDy^rMnOs, 
and was ascribed to an elliptical modulation of the spin 
cycloid [13, HO] ■ Figure [TT] shows the measured temper- 



ature dependence of y Sc{q)/Sb{q) for the 6c-plane spin 
cycloid in TbMnOs where q is the spiral propagation 
wave vector jToj . Importantly, the data for finite tem- 
peratures are extrapolated not to unity but to ^0.8 at 
T = 0. This may imply that the spin cycloid is ellip- 
tically modulated even in the ground state possibly due 
to quantum fluctuations. However, our calculation re- 
veals that the peaks of two spin correlation functions at 
T=0 are not equal even within our classical model [see 
insets of Figs. [T2^a) and fTSTb)]. in which the quantum 
fluctuations are not taken into account, indicating that 
the consideration of quantum fluctuations is not neces- 
sarily required to explain this experimental observation. 
In the following, we address this issue by investigating 
the detailed spin structures in the cycloidal spin states. 
We first show the calculated temperature dependence 

of y Sa{q)/Sb{q) for the a6-cycloidal case in Fig. [T2ja), 



and that of y Sc{q)/ Sb{q) for the 6c-cycloidal case in 
Fig. ini^b). They are calculated for J2= 0.80 meV and 
ac= 0.24 meV, and for 0.80 meV and ac= 0.38 
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FIG. 12: (Color online) (a) Temperature dependence of Sa{q), Sb{q) and y Sa{q)/ Sb{q) for the afe-cycloidal spin state with 

J2= 0.80 meV and a^— 0.24 meV. (b) Temperature dependence of Sb{q), Sc{q) and ^ S^q) / Sb{q) for the 6c-cycloidal spin 
state with J2= 0.80 meV and Uc— 0.38 meV. The data for finite temperatures obtained by the Monte-Carlo simulations are 
smoothly extrapolated to the data at T— obtained by numerically solving the Landau-Lifshitz-Gilbert equation as indicated 
by dashed lines. Insets of (a) and (b) show the spin-correlation functions in the k^-ky plane for fee = tt at T—Q. (c) Spin 
alignment (left panel) and spin directions (right panel) at r=0 for the afe-cycloidal spin state. Here the orthorhombic unit cell 
contains alternately stacked two different MnO planes with Zi—O and 0.5. (d) Those for the fec-cycloidal spin state. 



meV, respectively. Since the commensurate spin cycloid 
is convenient for discussion, we choose J2= 0.80 meV 
(i.e. Jab/ J2~ 1), for which the spiral propagation wave 
number becomes qm = 1/3. This value corresponds to 
an intermediate value of ~ 0.28 in TbMnOs and 
qm ~ 0.36 in DyMnOa. The data for finite tempera- 
tures are obtained by the Monte-Carlo simulations, while 
those for the ground states are obtained by numerically 
solving the Landau-Lifshitz-Gilbert equation. They are 
smoothly connected, and reproduce the experimentally 
observed temperature dependence and the non-unity val- 
ues at r=0 for both ab- and 6c-cycloidal cases. 

To clarify the origin of non-unity values of 



Sa{q)/Sb{q) and Sc{q) / Sb{q) at r=0, we investi- 
gate detailed spin structures in the ground state for both 
ah- and 6c-cycloidal cases. The results are shown in 
Figs. [TST c) and [T2r d). Noticeably, the rotation angles 
of cycloidal spins are not uniform, but suffer from con- 



siderable modulation. This rotation-angle modulation is 
due mainly to the single- ion spin anisotropics. The spins 
avoid to point along the hard magnetization axes, but 
prefer to direct along the easy magnetization axes. This 
is the very origin of the above inequivalent peak values of 
two kinds of spin correlation functions. We expect that 
higher harmonic peaks of the spin correlation functions 
should appear somewhere in the momentum space. In 
the a6-cycloidal case, we can see additional tiny peaks in 
the inset of Fig. [T^ a). 



E. Effects of single-ion anisotropies in the collinear 
spin states 

In the (nearly) collinear magnetic phases in i?Mn03 
like AFM(A) and sinusoidal phases, the Mn spins are 
aligned parallel to the h axis. The issue which inter ac- 
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FIG. 13: (Color online) (a) Local hard magnetization axes 
(solid arrows) and local easy magnetization axes (dashed 
arrows) at Mn A and Mn B sites due to the single-ion 
anisotropies, Tisia = Tifi^ + T^Sa- Owing to the tilting of 
MnOe octahedra, these axes deviate from the cubic x, y and 
z axes, (b) Polar representation of the 5=2 classical-vector 
spin. The angle (j>i is defined with respect to the orthorhom- 
bic a and b axes, [(c) and (d)] dependence of energies -E^^ai 
a^nd -Ecub at (c) Mn A and (d) Mn B for a spin sticking 
on the ab plane, [(e) and (f)] (j>i dependence of energies E^^^ji^, 
^SaB and -E^aA + ^SaB in the (e) presence and (f) absence 
of the MnOe tilting. Here -Ej'aA and -Ej'aB are the energies 
of H^a at Mn A and Mn B sites, respectively. For the MnOe 
tilting, we use the structural parameters of EuMnOa taken 
from Ref. 

tion or which magnetic anisotropy is responsible for this 
easy- axis spin anisotropy has not been clarified yet. In 
the Hamiltonian ([2]), there are some suspects, e.g. the 
single-ion anisotropy terms {Ti-^g, and H^^^), the DM- 
interaction term (TYdm), and the cubic- anisotropy term 
(Ticub)- Among them, the DM-interaction term turns 
out to be irrelevant. We confirm that the spins in these 
phases are still parallel to the b axis even without DM 
interaction. 

In the literature, Matsumoto discussed a combina- 
tion effect of cubic anisotropy Hcuh and one of the 
single-ion anisotropies Ti.^^ in the tilted MnOe octahe- 
dra [8^. Using the polar representation {Sai, Sbi, Sd) — 

{S sm9iC0S(j)i, S sm9ism(j)i, S cos9i) see Fig. [T^a). 

and considering that the hard magnetization along the c 



axis favors spins lying on the ab plane (i.e. 9i= tt/2), the 
term Ticub can be rewritten as 



'Wcub oc ''^^{S'^i + Sy^ + S'^^) 

i 

= 5'*^(3-cos40,)/4. (26) 



This expression implies that the cubic anisotropy favors 
spins pointing along the a or & axis since this term has 
energy minima at 0, 7r/2, tt, and 37r/2. Furthermore, 
tilting of the MnOe octahedra makes the local Q axes, 
which are local hard magnetization axes at every Mn 
site, inclined towards the a axis, resulting in a relatively 
hard magnetization along the a axis. Seemingly, this 
leads to a relatively easy magnetization along the b axis. 
The energy of Ticuh + 'Wfi'a polar representation 

with 9i — 'k/2 indeed gives energy minima when the spin 
points in the ±& direction (i.e. tt/2 and Stt/2). 

However, this scenario turns out to be wrong since the 
energies of Hcuh and (^cub and E^J are negligi- 
bly small as compared to the energy of {E^^) when 
the spins stick to the ab plane. Instead, the magnetic 
anisotropy is governed by the term Ti^^- Because of the 
term Tijfa: the £,i and rji axes become a hard magnetiza- 
tion axis alternately in the ab plane. In Figs. [T^T c) and 
[TSf d). the ipi dependence of E^^ is shown together with 
those of -Ecub and E^^ for Mn A and Mn B sites, respec- 
tively. Here Mn A and Mn B are two different Mn sites in 

the same Mn-0 plane see Fig.[lja). Noticeably, the 

energies -Ecub and E^^ are much smaller than the energy 
E^^. However, the energy E^^{(j)i) has energy minima 
not along the b axis {4ii = tt/2), but nearly along the x 
and y axes ((/)i=7r/4 and 37r/4), respectively. This may 
seem incompatible with the easy magnetization along the 
b axis. However, since the in-plane ferromagnetic ex- 
changes Jab predominantly force the spins at Mn A and 
Mn B sites to be parallel, a sum of these two energies, i.e. 
E|, at Mn A and E|, at Mn B (Ef^^ and E^^) gov- 
erns the magnetic anisotropy in the coUinear spin phases. 
We show the (f>i dependence of energies E^^j^, -E'ifaB and 
^^sfaA+B=^^sfaA+^^laB in Fig. [He). This figure actually 
shows that the energy minimum of E^^j^^Q{(f>i) is located 
at the b axis (</>,; =7r/2). 

In fact, this is a consequence of the MnOg tilting. 
Without tilting, the energy sum E^g^^^g((/)i) is always 

zero irrespective of (pi see Fig. [I^f). Although both 

^SaAi'Pi) and -Ejfg^g((/)i) have strong (pi dependence, they 
perfectly cancel out. The above discussion is valid also for 
Mn C and Mn D sites in another Mn-0 plane because of 
the presence of a mirror plane between two Mn-0 planes. 
We conclude that the combination effect of the single-ion 
anisotropies Hsia = Wf?a + ^^a and the GdFeOs-type 
distortion is an origin of the easy-axis spin anisotropy in 
the AFM(A) and sinusoidal phases. 
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FIG. 14: T-E phase diagram obtained by the Monte-Carlo 
analysis of the Hamiltonian ([2|. Here the parameter _B is a 
variable which expresses the strength of single-ion anisotropy 
Tifia, and J2 is fixed at 0.80 meV. Values of other model pa- 
rameters in ([2]) are fixed at the values listed in Table IIIII (see 
Sec. |llj. The temperature range over which the sinusoidal 
coUinear spin phase emerges increases as E is increased. 



F. Sinusoidal collinear spin phase 
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FIG. 15: T-K phase diagram obtained by the Monte-Carlo 
calculation for the case without GdFeOa-type distortion (see 
text). The term Ti.K given in Ea. (|27|l is added to the Hamilto- 
nian (121), which generates the easy-axis spin anisotropy along 
the b axis. Here, the parameter if is a variable which ex- 
presses the strength of easy-axis spin anisotropy. Values of 
the model parameters except for K and J2=0.80 meV are 
fixed at the values listed in Table [HI] (see Sec.|IT]|. The tem- 
perature range of the sinusoidal collinear spin phase increases 
as K is increased. 



In this section, we discuss the nature and origin of the 
sinusoidal coUinear spin phase in the intermediate tem- 
perature regime. In this phase, the spins are ahgned 
along the b axis with sinusoidally modulated amplitudes. 
We find that for its emergence, the single-ion anisotropy 
term 7^^^^ is essentially important. In Fig. [HI we display 
evolution of each magnetic phase as a function of the 
anisotropy parameter E. This figure shows that the tem- 
perature range over which the sinusoidal collinear phase 
emerges increases as E is increased. 

As we have discussed in the previous section, when the 
GdFeOs-type distortion is present, the term Hf^^ gener- 
ates an easy-axis spin anisotropy along the b axis in the 
collinear spin phase. In fact, this easy-axis anisotropy 
is decisively important for the stability of the sinu- 
soidal collinear phase. To see the importance of easy- 
axis anisotropy, we have performed calculations for the 
cubic-symmetry case, i.e. the case without GdFeOs-type 
distortion. We set a—b, a;oi=0, ?/Oi=0.5, a;o2=0-75, 
2/02=0-25, and 2:02=0 for the local-axis vectors ^i, rji 
and in Eqs. ([T])-®. As a result, the local axes ^i, rji 
and Q attached to the i-th MnOg are equivalent to the 
crystallographic x, y and z axes. In this case, the term 
Tisia does not generate an easy-axis anisotropy, and the 
sinusoidal collinear phase indeed disappears. Instead, we 
further add a term of the easy-axis anisotropy Hk by 
hand to the Hamiltonian which is given by 

Hk^-KJ^SI, (27) 

i 

with K > 0. We obtain similar evolution of the sinusoidal 

collinear phase as a function of the parameter K see 

Fig. [121 Note that when K=0, the sinusoidal collinear 



phase vanishes, indicating that the GdFeOa-type distor- 
tion or the easy-axis anisotropy is indispensable for its 
emergence. We indeed notice that in the experimental 
phase diagram of Fig. [IJb), the temperature range of 
the sinusoidal collinear phase is larger in a material with 
stronger GdFeOs-type distortion, or a smaller R ion. 

The above facts can be understood easily as previously 
discussed in Ref. [l05| |. As temperature is lowered, the 
first ordered state has spins aligned along the easy axis, 
i.e. the b axis. As temperature is further lowered, the 
spin length grows. Then, the system develops a long- 
range transverse sinusoidal order to satisfy the compe- 
tition between Jab and J2. This requires to enter the 
cycloidal spin phase. The temperature range of the si- 
nusoidal collinear phase increases as the anisotropy K is 
increased because the longer spin length is required to 
overcome the larger K, which in turn requires going to 
lower temperature. 



IV. CONCLUSION AND DISCUSSION 

In summary, we have studied the magnetoelectric 
phase diagrams of the perovskite manganites /JMnOa by 
constructing a microscopic spin model to describe the Mn 
3d-spin systems. We have analyzed this model by using 
the Monte-Carlo method for the thermodynamic prop- 
erties and by numerically solving the Landau-Lifshitz- 
Gilbert equation for the ground-state properties. Consid- 
ering that the GdFeOs-type lattice distortion enhances 
the second-neighbor antiferromagnetic exchanges J2, we 
have studied the T-J2 phase diagrams, and have obtained 
diagrams in good agreement with the experimental ones 



including two kinds of multiferroic phases, i.e. the ab- 
cycloidal spin phase with Pa and the 5c-cycloidal spin 
phase with Pc- 

We have discussed a mechanism of the electric polar- 
ization flop. The afo-cycloidal spin state is stabilized by 
the single-ion anisotropy H^^^ and the DM interaction 
Tigfj^ with vectors on the in-plane bonds. Here the single- 
ion anisotropy Ti^j^ makes magnetization along the c axis 
hard. On the other hand, the 6c-cycloidal spin state 
is stabilized by the DM interaction Hfy^ with vectors 
on the out-of-plane bonds. As the spiral rotation an- 
gle increases with increasing J2, the energy gain due to 
'^DM is reduced in the a5-cycloidal spin state. In this 
way, the increasing GdFeOa-type distortion destabilizes 
the a6-cycloidal spin state. This leads to the cycloidal- 
plane flop from ab to be, and consequently to the electric- 
polarization flop from Pa to Pc- 

We have demonstrated that the regime of 6c-cycloidal 
spin phase with Pc increases with increasing ac, and have 
revealed that the phase diagram of Gdi_2;Tba;Mn03 with 
a large Pc regime is reproduced for a rather large value of 
ac— 0.38 meV, while that of Eui-ajYj^MnOs with a small 
Pc regime is reproduced for a smaller value of ac— 0.30 
meV. So far existence of the two different ferroelectric 
phases with Pa and Pc in KMnOs compounds have often 
been connected with the role of /-electron moments on 
the rare-earth ions. However our result indicates that the 
rare-earth magnetism plays a relatively minor role. 

By investigating the detailed structures of the cycloidal 
spin states, we have found that the spiral rotation an- 
gles are not uniformly the same but are significantly 
distributed due to the single-ion anisotropics. The ex- 
perimentally observed direction-dependent magnitude of 
spin-correlation peak is ascribed to this rotation-angle 
distribution instead of the elliptical modulation of the 
spin cycloid as claimed so far. 

By examining the effect of single-ion anisotropics, we 
have found that because of the tilting of MnOg octa- 
hedra, the single-ion anisotropics, Hsia = Wfi'a + ^^a' 
energetically favor spins pointing along the b axis in the 
AFM(A) and sinusoidal coUinear phases. 

We also found the sinusoidal collinear spin phase is 
stabilized by the easy-axis spin anisotropy along the b 
axis generated by the single-ion anisotropy term Hsia. in 
the distorted lattice structure of GdFeOa type. 

The microscopic model proposed here must be helpful 
for studying origins and mechanisms of several intrigu- 
ing magnetoelectric phenomena discovered in the present 
manganite compounds, e.g., origins of t he el ectrically 
activated spin excitation (electromagnon) [l06| |. mecha- 
nism s of the magn etic-field-induced electric polarization 
fiop |l07l . llOSl . llOOll. and origins of the giant magnetoca- 
pacitance effect [411 l42j . Discussion and analysis of the 
spin-wave di spers i ons obtained in the neutron-scattering 
experiments [llOl 111 j | are also an issue of interest. 

Now, we would like to compare our work with the re- 
cent theoretical work done by Dagotto and coworkers, 
which is based on the two-orbital double-exchange model 
with additional terms [III [s^ . In Ref. [Ill ' they added 
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FIG. 16: (Color online) T hird-nelRhbor spin exchange J3 con- 
sidered in Refs. [ll5l . lll^ . 



the DM-interaction term with DM vectors coupling to 
the lattice, and reproduced the cycloidal spin order. On 
the other hand, in Ref. [3^, they considered very weak 
second-neighbor spin exchanges to reproduce a succes- 
sive emergence of the AFM(A), cycloidal, and AFM(E) 
phases, in the ground state. However, in those studies, 
they failed to reproduce the cycloidal-plane flop and the 
sinusoidal collinear phase. In addition, within their the- 
ory, it seems hard to explain the monotonical increase 
of the propagation wave number g^" (or the spiral ro- 
tation angle 9s) with increasing GdFeOs-type distortion 
observed in experiments [3|40^[83]- On the other hand, 
we consider that the second-neighbor exchanges J2 play 
a predominant role in stabilizing the cycloidal spin state, 
and the ferroelectric polarization subsequently emerges 
through the inverse DM mechanism. Our model in which 
the strength of J2 scales with the orthorhombic lattice 
distortion naturally explains the observed monotonical 
increase of g^". Moreover, the spiral rotation angle de- 
termined by the ratio Jab/J2 turns out to control the 
competition between the ab- and 6c-cycloidal spin states, 
resulting in the electric-polarization flop. This indicates 
that consideration of the second-neighbor exchanges J2 
is essentially important to describe the magnetoelectric 
coupling in the manganite system. Importance of the 
further exchanges is peculiar in the present manganites. 
This is owing to the weak nearest-neighbor exchanges re- 
sulting from the cancellation of ferromagnetic and anti- 
ferromagnetic contributions from the Cg- and i2g-orbital 
sectors due to the ^29^3 electron confi guration. Indeed 
the values of Jab ^1.6 meV in LaMnOs Ill2llll3 ll is much 
smaller than Jab— 13.3 meV in LaTiO-? |ll4l |. which is a 
S—1/2 G-type antiferromagnet with electron config- 
uration. 

It is also worthy mentioning that on the basis of the 
LDA-I-U calculation, Solovyev proposed that the third- 
neighbor spin exchange J3 (see Fig. is increased 
by the GdFeOa-type distortion, and this exchange also 
plays an impor tant role for stability of the spiral spin 
order |ll5l . Illdj . In this case, the propagation wave vec- 
tor of the magnetic spiral is determined by the frustra- 
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tion among three exchanges Jab^ Ji and J3, and its wave 
number is given by, 



1 _l r '^ah n 



(28) 



To check whether the consideration of J3 affects the re- 
sults presented in this paper, we have examined the ef- 
fects of J3. We have found that once the wave number 
qm is determined, the same resuhs are obtained. 
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APPENDIX A 

In this Appendix, we calculate the values of superex- 
change parameters, J ah and Jc, for /JMnOa by fol- 
lowing the formulation given by Gontchar and cowork- 
ers 82, 83]. 

In 7?Mn03, the ground-state orbital state of each 
Mn'^+ ion is dominantly determined by local coordina- 
tions of surrounding oxygens with the Jahn- Teller distor- 
tion because of the strong electron-lattice coupling. The 
wave functions of occupied and hole orbitals at the n-th 
Mn^+ site, and ^'2n, can be written as linear com- 
binations of the twofold eg-orbital wave functions [ll7l | ; 



*i„ = sin— (^„e -f cos— V3„e, 
= COS— (p„e - sin— V3„£, 



with 



(Al) 
(A2) 

(A3) 
(A4) 



where rjn and Qn are coordinates with respect to the 
local axes attached to each MnOg octahedron. 

The values of orbital angle 0„ can be derived from ex- 
perimental data, using parameters of lattice distortions; 



Qe 



smt 



Qe 



(A5) 
(A6) 



with Qon and being symmetrized Jahn- Teller distor- 
tions on each MnOg octahedron of 3z^ — and — 
types, respectively. 



In the orthorhombic i?Mn03 crystal described by the 
Pbnm space group, the Jahn- Teller distortions Qgn and 
Qen are given by 



Qe = 



Qe = V2ivxa + Vyb), 



with 



-^-^02, 



4 -yo.. 



(A7) 
(A8) 

(A9) 
(AlO) 



Here, Vx and Vy are shifts of the in-plane oxygen envi- 
ronment, XO2 and 2/O2 are coordination parameters for 
the in-plane oxygen ions, and a, b and c are the lattice 

parameters in Pbnm axes see Table. IIVI 

The Jahn- Teller distortion in i?Mn03 causes a C-type 
orbital ordering, for which the orbital angles satisfy the 
following relation: 

(/-i = 9!'2 = -h = -(l^i = 0- (AH) 

Within this framework, the superexchange parameters 
depend on the lattice distortions and the orbital structure 



Jab = '°?o [1 - « COS + /3(cos^ + ^)IA12) 
J, = -I^^^^ll + 2a COS cl)~ 13 cos^^], (A13) 



where ipab and ipc are means for the in-plane and out-of- 
plane Mn-O-Mn bond angles, respectively, and Vab and Tc 
are means for the in-plane and out-of-plane Mn-0 bond 
lengths, respectively. For the values of parameters ap- 
peared in these formulas, we use Jn=1.456 eVA, a=1.0 
and (3=4.5 according to Refs. fs^. [ssj. 

The values of Jab and Jc calculated using struc- 
tural data in Refs. [H, [H IM 113, Ull are plotted in 
Fig. [T7| as functions of the ionic i?-site radius. Here, 
the values for LaMn03 were obtained from the spin- 
wave dispersion measured in the neutron-scattering ex- 
periments as Jab= 1-66 meV and Jc= 1.16 meV in 
Ref. Ill2ll . and Jab— 1-67 meV and Jc= 1.21 meV in 
Ref. [1135. The values for PrMn03 were also obtained in 
the neutron-scattering experim ent a s Jab= 1.12-1.19 meV 
and Jc= 1.2-1.29 meV in Ref. [ull- The calculated val- 
ues for LaMn03 and PrMn03 are in good agreement with 
these experimental values. We note that both Jab and Jc 
are nearly constant in the region of relatively small-sized 
R ions, e.g. i?=Eu, Tb, Dy and Y, in/near which the 
multiferroic phases emerge. 



APPENDIX B 

In this Appendix, we estimate the values of parameter 
J2 for i?Mn03. The spiral spin orders in i?Mn03 have 
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TABLE IV: Structural parameters for .RMnOs. Here r^j and rj;, denote lengths of the longer and shorter in-plane Mn-O bonds, 
respectively. References from which the data are taken are shown in the second column. 



RMnOs 


Ref. 


IR(A) 


a{A) 


b{A) 


c{A) 




xo^ 




yo2 




Z02 


(j}{d 


sg) 


fc{deg) 


V5a6(deg) 


rc{A) 






LaMnOs 


([84]) 


1 


216 


5.5367 


5.7473 


7.6929 





7256 





3066 





0384 


107 


06 


155.48 


155.11 


1.9680 


2 


178 


1 


907 


LaMnOa 


([87]) 


1 


216 


5.5405 


5.7458 


7.6998 





7261 





3069 





0388 


106 


43 


155.28 


155.04 


1.9706 


2 


184 


1 


904 


PrMnOa 


([85]) 


1 


179 


5.4491 


5.8129 


7.5856 





7151 





3174 





0430 


110 


93 


152.36 


150.51 


1.9530 


2 


210 


1 


909 


PrMnOa 


([83) 


1 


179 


5.4500 


5.8295 


7.5805 





7151 





3179 





0433 


111 


41 


151.59 


150.36 


1.9549 


2 


217 


1 


910 


NdMnOa 


([85]) 


1 


163 


5.4170 


5.8317 


7.5546 





7141 





3188 





0450 


111 


52 


150.86 


149.63 


1.9514 


2 


218 


1 


905 


NdMnOa 


(M) 


1 


163 


5.416 


5.849 


7.543 





7132 





3204 





0450 


112 


09 


150.3 


149.2 


1.951 


2 


227 


1 


905 


NdMnOa 


([87]) 


1 


163 


5.4168 


5.8518 


7.5479 





7124 





3199 





0447 


112 


67 


150.38 


149.27 


1.9520 


2 


223 


1 


911 


SmMnOa 


([86]) 


1 


132 


5.362 


5.862 


7.477 





7076 





3241 





0485 


114 


37 


147.6 


147.0 


1.947 


2 


232 


1 


910 


EuMnOa 


([86]) 


1 


120 


5.340 


5.866 


7.448 





7065 





3254 





0487 


114 


82 


146.54 


146.47 


1.944 


2 


234 


1 


907 


EuMnOa 


([87]) 


1 


120 


5.3437 


5.8361 


7.4619 





7055 





3247 





0485 


114 


55 


147.35 


146.45 


1.9438 


2 


220 


1 


912 


GdMnOa 


([86]) 


1 


107 


5.318 


5.866 


7.431 





7057 





3246 





0508 


115 


74 


145.6 


146.0 


1.944 


2 


229 


1 


915 


TbMnOa 


([85]) 


1 


095 


5.2931 


5.8384 


7.4025 





7039 





3262 





0510 


115 


47 


145.06 


145.36 


1.9401 


2 


221 


1 


905 


DyMnOa 


(^5J) 


1 


083 


5.2785 


5.8337 


7.3778 





7028 





3276 





0521 


115 


69 


143.23 


144.70 


1.9437 


2 


224 


1 


903 


DyMnOa 




1 


083 


5.2802 


5.8448 


7.3789 





7013 





3249 





0518 


118 


86 


145.37 


145.03 


1.9322 


2 


210 


1 


919 


YMnOa 


/ To rl \ 

([85]) 


1 


075 


5.2418 


5.8029 


7.3643 





7005 





3266 





0520 


116 


17 


143.51 


144.55 


1.9385 


2 


200 


1 


904 


HoMnOa 


([85]) 


1 


072 


5.2572 


5.8354 


7.3606 





7013 





3281 





0534 


116 


28 


142.47 


144.08 


1.9435 


2 


222 


1 


905 


ErMnOa 


([85]) 


1 


062 


5.2262 


5.7932 


7.3486 





7003 





3266 





0542 


116 


18 


142.82 


143.92 


1.9382 


2 


199 


1 


903 


ErMnOa 


([88]) 


1 


062 


5.2395 


5.8223 


7.3357 





6997 





3343 





0563 


113 


59 


141.59 


142.05 


1.942 


2 


248 


1 


891 


TmMnOa 


([88]) 


1 


052 


5.2277 


5.8085 


7.3175 





6984 





3372 





0571 


112 


75 


141.15 


141.17 


1.940 


2 


255 


1 


886 


YbMnOa 


([88]) 


1 


042 


5.2163 


5.7991 


7.2992 





6979 





3394 





0580 


112 


07 


140.27 


140.51 


1.940 


2 


263 


1 


879 


LuMnOa 


([88]) 


1 


032 


5.1972 


5.7868 


7.2959 





6989 





3415 





0575 


110 


14 


139.56 


140.36 


1.944 


2 


269 


1 


862 



Structure data from 
-^■o- Alonso et al.* 

Mori etal. 
-^-n- Dabrowski et al. 

Tachibana et al. 
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17: (Color online) Calculated superexchange parame- 
Jab and Jc, plotted as functions of the ionic _R-site ra- 
Structural data used in the calculation are taken from 
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FIG. 18: (Color online) Estimated values of J2 in .RMnOa and 
several solid solutions plotted as functions of the ionic i?-site 
radius. Inset describes the J1-J2 classical Heisenberg model 
on the anisotropic triangular lattice and the spiral spin order 
with uniform rotation angles of 9. 



turned out to be realized due to the competition between 
the nearest-neighbor ferromagnetic couphng Jab and the 
second-neighbor antiferromagnetic coupling J2 in the ab 
plane. Along the c axis, the spins are forced to be aligned 
antiferromagnetically because of the strong antiferromag- 
netic coupling Jc- Thus the propagation wave number 
and the rotation angle of the magnetic spiral are domi- 
nantly determined by the ratio J2/ Jab although magnetic 



anisotropics and DM interactions should modulate them 
to some extent. Under this circumstance, we can roughly 
estimate the values of J2 with the aid of the simple J1-J2 
classical Heisenberg model on the anisotropic triangular 
lattice (see inset of Fig. [TS]) . The Hamiltonian is given 
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by 



<i,j> 



Here, Ji and J2 are the ferromagnetic exchange for the 
nearest-neighbor bonds and the antiferromagnetic ex- 
change for the second-neighbor (diagonal) bonds, respec- 
tively. We introduce the third term with K > to deter- 
mine the basal spiral plane by making the magnetization 
along the z axis hard. 

Within this model, the energy for spiral spin state per 
site is given as a function of the spiral rotation angle 9; 



E{0)/S^ = -2 Ji cos 9 + J2 cos(26'). 



(B2) 



From a saddle-point equation, dE{6)/dd = 0, we obtain 
a relation cos^s = Ji/(2J2) where 9s is the rotation an- 
gle of spiral spin state with a minimum energy. This 
formula implies that a spiral spin state emerges when 
J2/J1 > 0.5, while when J2/J1 < 0.5 a ferromagnetic 
state is stabilized. 

We estimate the values of J2 in TbMnOs, DyMnOa, 
YMnOa and several solid solutions with the use of exper- 
imentally measured spiral wave numbers g^" \2Mt, llO, l80|] 
and the values of Ji {Jab) calculated in Appendix A (see 
Table W} . In Fig. [181 we plot them as a function of the 
ionic i?-site radius, which exhibits a monotonical increase 
with decreasing i?-site radius. We can deduce the values 
even for RMnO^ with i?=Eu and Gd by extrapolating 
the data, although their ground-state magnetic structure 
is not spiral but AFM(A). 

When g^" is small, the energy difference between the 
spiral spin state and the AFM(A) state would be so small 
that the DM interaction and the single-ion anisotropy can 
affect their energetics significantly. Both of these two fa- 
vor the AFM(A) state. The collinearly aligned spins par- 
allel to the b axis in the AFM(A) state can benefit from 
both the easy magnetization b axis due to the single-ion 
anisotropy and the DM vectors on the out-of-plane bonds 
with a large a component ac- Therefore, the estima- 
tion of J2 based on the pure J1-J2 model performed here 
would involve ambiguities particularly near the AFM(A)- 
cycloidal phase boundary or in TbMnOa. The actual 
value of J2 in TbMnOs is expected to be slightly larger 
than the above-calculated value of J2 ~ 0.62 meV. 



APPENDIX C 

In this Appendix, we calculate the values of single-ion 
anisotropy parameters D and E in i?Mn03 with the use 
of formulas proposed by Matsumoto [soj . 

The parameters D and E can be calculated within the 
second-order perturbation theory in terms of the spin- 
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FIG. 19: (Color online) Calculated single-ion anisotropy pa- 
rameters, D and E plotted as functions of the ionic _R-site 
radius. Structural data used in the calculation are taken from 
Refs. [sj. Issl. Is^. [stI. IssI] . Optical absorption data are taken 
from different references between (a) and (b), i.e. Ref. [89l | 
for (a) and Ref. ;90] for (b). 



orbit interaction. They are given by 



D 



-3( 



A2 



4A^ 



AS 3A 



) cos 4>„ 



E = : sm (b„ , 

AE ^ ' 



(CI) 
(C2) 



where A is the coupling constant of the spin-orbit inter- 
action, and 0„ is the orbital angle which represents the 
orbital wave functions at the n-th Mn'^^ site as defined 
in Appendix A. Here, AE and A are energy levels for 
the two excited states of Mn^+ ion, ^Tig (de^) and ^T2g 
(de^c?7^), measured from the ground-state ^ 
level, and satisfy the following relation |l2 



Eg {de^dj) 



AE^6B + 5C-A, 



(C3) 



where B and C are the Racah parameters. Using the val- 
ues for free Mn'^+ ion, J5o=120 meV and Co=552 meV, we 
can obtain the orbital reduction factor x from Eq. (|C3p 
by setting B = xBq and C = xCq. Then, the value of A 
in the iJMnOa compounds can be obtained from A = xXq 
with Ao=10.5 meV being the couphng constant for free 
Mn^+ ion. 

The values of AE and A for Mn^+ ion surrounded 
by octahedrally coordinated oxygens are experimentally 
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TABLE V: On the basis of the simple two-dimensional J1-J2 classical Heisenberg model, the values of J2 in TbMnOa, DyMnOs, 
YMnOa and several solid solutions are estimated from the experimental data of the spiral wave numbers gj^" and theoretically 
calculated values of Jab- Here, IR denotes the ionic radius of the R ion, and 6 is the spiral rotation angle calculated from gj^" 
as 9 — 180° X Qm". The references from which the data of are taken are presented in the last column. The effective _R-site 
radii in the solid solutions are deduced by interpolation. 



iJMnOs 


IR(A) 2q^^ 


e (deg 


) Ji (meV) 7 = 


(2 cos 6) 


J2 =7J'i ( 


meV) Ref. 


TbMnOa 


1.095 0.56 


50.4 


0.79 


0.784 


0.62 


[40] 


DyMnOs 


1.083 0.78 


70.2 


0.76 


1.47 


1.12 


[40] 


YMnOa 


1.075 0.87 


78.3 


0.77 


2.466 


1.90 


[119] 


Euo.6Yo.4Mn03 


1.102 0.58 


52.2 


0.817 


0.815 


0.67 


[18] 


Tbo.32Dyo.68Mn03 


1.087 0.678 


61.0 


0.768 


1.032 


0.79 


[80] 


Tbo.4iDyo.59Mn03 


1.088 0.663 


59.6 


0.786 


0.990 


0.78 


[80] 


Tbo.5oDyo.5oMn03 


1.089 0.656 


59.0 


0.773 


0.972 


0.75 


[80] 



TABLE VI: Energy levels AS and A for Mn^+ ion surrounded 
by octahedrally coordinated oxygens obtained in the optical 
absorption measurements from Refs. [90I]. Calculated or- 
bital reduction factor x and coupling constant of the spin- 
orbit interaction A are also listed. 



A£; (meV) A (meV) 


X 


A (meV) Ref 


353.4 2530.8 
349.7 2083.2 


0.83 
0.71 


8.75 [89] 
7.38 [90] 



obtained in the optical absorption spectra [H, [o^. In 
Table. IVI) we list the values of /S.E and A taken from 
two different references as well as the values of x and 
A. Usin g th ese parameter values and structural data 
in Refs. [H, HI, HE [13, HI], we calculate D and E 



for i?Mn03 with various R ions, which are plotted in 
Figs. fTOT a) and fTOT b) as functions of the ionic i?-site ra- 
dius. Note that the difference of A£' and A between two 
references causes a small difference of calculated values 
of D and E. In our model calculation, we use the val- 
ues obtained from the optical data in Ref. [93|, because 
with these values, the calculation reproduces the experi- 
mental phase diagrams with a reasonable set of the DM 
parameters, which is consistent with the band calcula- 
tion 94] and the ESR experiments [9^, [o^] as discussed 
in Sec. Ill Al On the other hand, with the values ob- 
tained from the optical data in Ref. [s^, we need to as- 
sume rather large DM vectors on the out-of-plane bonds 
to reproduce the experimental phase diagrams. 
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